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Abstract 

We report the results of monitoring the four images of the gravitational lens 
MG 0414+0534 at 15 GHz. In 35 VLA maps spanning 180 days, we measure 
root-mean-square variations in the image light curves of ^3.5% mostly due to 
variations in the flux density calibration. The flux ratios, which are indepen- 
dent of flux density calibration variations, show root-mean-square variability 
of 1-3%. Extensive simulations of the data analysis process show that the ob- 
served variations in the flux ratios are likely to be due entirely to errors in the 
deconvolution process. It is possible that some of the observed variation is due 
to the source; however, the signal-to-noise ratio is too small to make a time 
delay determination using a data set of this size. 
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1 Introduction 



Well before the first gravitational lens was discovered, Refsdal ( |1964a| ; |1964b ) pointed 
out that, given a model of the lensing potential and the redshift of both the source 
and the lens, a measurement of the difference in optical path length between two 
images in a gravitational lens could be used to determine the Hubble parameter. 



More recently, Narayan ( 1991 ) has shown that a time delay measurement combined 
with a model of the lensing potential provides a measure of the angular diameter 
distance to the lens that is independent of the redshift of the source and cosmolog- 
ical assumptions other than local isotropy and homogeneity transverse to the line 
of sight. Thus, lens models combined with time delay measurements and lens red- 
shifts can be used to measure the relationship between angular diameter distance 
and redshift. The measurement of the angular diameter distance-redshift relation 
at a variety of lens redshifts will allow the comparison of the measured relation to 
that predicted by cosmological models. Large-scale structure may contribute to the 



deflection and time delays in strong lensing ( Seljak 1994 ; Bar-Kana 1996 ) raising 
the interesting possibility of detecting the effect of large-scale structure in a sample 
of lenses. 

For a time delay measurement we have selected the gravitational lens MG 0414+0534, 
a source that was culled from the MIT-Greenbank survey ( |Bennett et al. 1986 ) in 
a gravitational lens search ( Hewitt et al. 1988[) and followed up with optical and 
radio observations. The object consists of four bright, highly reddened images with 
very similar radio spectra and optical colors. These characteristics, along with 
the morphology and relative brightnesses of the images, were taken as early evi- 
dence that the object is indeed a gravitational lens ( [Hewitt et al. 1992j ), Further 
study has supported that finding. Katz and Hewitt ( |1993| ) confirmed that the close 
double (components Al and A2) is indeed a double as predicted by lens models. 
Optical (I-band) observations at the Michigan-Dartmouth-MIT observatory at Kitt 
Peak have detected the lensing galaxy at nearly the expected position ( [Schechter fc 
Moore 1993] ). Spectroscopy by Angonin-Willaime et al. ( 1994| ) shows that the spec- 
tra of A1+A2 and B are similar. Infrared spectroscopy by Lawrence et al. ( |1995| ) 
suggests that the source is a fairly typical quasar that is highly reddened by dust in 
the lensing galaxy (but see Vanderriest et al. (1995) and Annis et q/.( |1993 ) for other 
views). A 15 GHz radio image of MG0414 is displayed in Figure |l[ 

MG0414 is well suited to relating time delay measurements to angular diameter 
distances since the morphology of the system argues for a relatively simple lensing 
potential (particularly compared to QSO B0957+561). Furthermore, HST, VLA, 
and VLBI studies of the object are underway which have revealed structure that may 
be used to place constraints on models of the lensing potential (|Falco et al. 1997 ; 
Trotter Wilson 1997| ; |Katz et al. 1997j ; |Ellithorpe 19951 ; |Patnaik & Porcas 1995|)~ 
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Figure 1: A 15 GHz (A-array) map of MG0414 produced in the course of the 
present work (observation date 15 Dec 1992). The contours correspond to flux 
densities of 0.8, 1.6, 3.2, 6.4, 12.8, 25.6, 51.2, and 102.4 mJy. The circle in the 
lower right corner illustrates the size of the synthesized beam for this observation 
(.15"x.l4"). 



1.1 Lens Models 

Kochanek (1991) has fit five simple models of the lensing potential in MG0414 using 
only the positions of the images as constraints for reasons of computational speed 
and lack of reliable data on the flux ratios. The models fall into two classes: those 
in which image B is at a minimum of the time delay surface (models 2, 3, and 4) 
and those in which it is at a saddle point (models 1 and 5). In order to illustrate 
the difference between the two classes of models we show Kochanek's model 1 (Sin- 
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gular Isothermal Sphere + Internal Quadrupole) and model 3 (Singular Isothermal 
Sphere + External Quadrupole). Since Kochanek does not publish the fitted source 
positions, we estimate the source position by computing the source position implied 
by each observed image. Images Al, A2, and B give consistent estimates of the 
source position and we take their average as our estimate (0.053",— 0.050" from the 
lens center for model 1 and —0.13", 0.046" for model 3). Figure || shows the lensing 
potentials for both models (from Kochanek's fits) and our estimated source position. 
Even though the position angles of the potentials differ by approximately 90°, the 




Figure 2: The gravitational potential for Kochanek's model 1 (left) and model 3 
(right). The contours are equipotentials which are equivalent (in the weak field 
limit) to lines of constant index of refraction. The diamond signifies our estimated 
source position. The central detail in model 1 has been suppressed in order to make 
the source position visible. 

radial deflections and therefore the Einstein rings generated by these models are 
nearly aligned. 

Figure || displays the time delay surfaces corresponding to the two models. By 
the Fermat principle, images are formed at the extrema (marked) of the time delay 
surface. Both models can reproduce the positions of the images, but they do so 
very differently: saddle points and minima are exchanged in the two time delay 
surfaces. The B:C flux ratio predicted by Kochanek's models 1 and 5 (.58 and .51 
respectively) is inconsistent with our observations where the measured flux ratio is 
2.53. Models 2, 3, and 4 come closer to the measured B:C flux ratio (2.2-3.2) but 
fail to correctly predict the Al:A2 flux ratio. The very large magnification gradient 
in the neighborhood of Al and A2 should make the flux ratio of A1:A2 difficult to 
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Figure 3: The time delay surface for Kochanek's model l(left) and model 3(right) 
using the assumed source positions shown in Figure ||. Extrema (with the exception 
of the central maximum) are labelled "M" for minimum and "S" for saddle point. 

model. 

The models that are consistent with the measured flux ratios predict that flux 
variations should be observed first in image B, following next in Al & A2, and 
finally in image C. The model of the system by Hewitt et al. (1992) predicts (for 



Hq = 80 kms _1 Mpc _1 , qo = 0.5, empty beam) that images B and A2 lead Al by 
12 and 0.3 days respectively and that image C lags Al by 19 days. It should be 
noted that the model upon which these estimates are based should be taken to be 
illustrative only since it assumes a lens redshift (zi = 0.47) and does not exactly 
reproduce the observed geometry and flux ratios of the system. 

2 Observations 

We observed MG0414 at 15 GHz with the Very Large Array^ (VLA) starting on 
1992 November 8 in its most extended (A-array) configuration, continuing through 
the second most extended (B-array) configuration, and ending on 1993 May 3. The 
central frequencies for the two observing bands were the standard VLA U-band val- 
ues of 14.9649 and 14.9149 GHz, each with 50 MHz bandwidth. Each band was ob- 
served in both senses of circular polarization bringing the total observed bandwidth 
to 100 MHz in each polarization. There were a total of 63 individual observations 



3 The VLA is part of the National Radio Astronomy Observatory, which is operated by Associated 
Universities, Inc., under cooperative agreement with the National Science Foundation. 
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each of 0.5 to 1.5 hours duration. The total time spanned by the observing program 
and the sampling interval of the observations were chosen so that for a reasonable 
range of the (unknown) lens redshift, at least one delay would be detectable. This 
program was scheduled into small interstices in the VLA schedule resulting in obser- 
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Table 1: Journal of Observations 

vations occurring over a wide range of hour angle and with no particular sampling 
pattern (see Table |l]). We have found that, in practice, an irregular sampling pattern 
is desirable because it minimizes windowing effects in the time delay analysis. 

Some of the observations took place when the source was at very low elevation. 
These data are problematic since the projection of the baselines of the array and 
the resulting large synthesized beam makes deconvolution difficult, and because the 
atmospheric contribution to phase errors and flux calibration errors is much larger. 

3 Data Reduction 

The data were reduced using the National Radio Astronomy Observatory's Astro- 
nomical Image Processing System (AIPS). After the initial excision of raw data 
that were clearly corrupted by interference or hardware problems, all data sets were 
processed by AIPS with a "run file" that applied the same processing to each ob- 
servation. Complex antenna gains were determined by observations of J0423— 013. 
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Refined estimates of the antenna phases were obtained by self-calibration (Cornwcll 
fe Fomalont 1989Q . Observations of 3C 48 were used to set the flux density scale but 
in order to use as many of the available baselines as possible, we used a model of the 
source based on four minutes of A- array observation kindly provided by C. Katz. 

We used the CLEAN (Hogbom 1974) algorithm implemented within the AIPS 
task MX ( |Clark 1980 ; Schwab 1984) to deconvolve each image. CLEAN is an iter- 
ative procedure that, for each iteration, selects the brightest pixel in the map and 
subtracts a fraction of it multiplied by the dirty beam. As a result, when two point 
sources are very close together, components tend to fall spuriously between them. 
This effect becomes particularly troublesome for the B-array data in which the Al- 
A2 doublet is barely resolved. Our initial attempts to clean these data resulted in 
a single bright source located between the two images with flux extending from the 
center towards the true image positions. When we instead used the knowledge of 
the relative positions of the components as measured from our A-array maps to con- 
strain CLEAN so that it only attempted to place clean components in 0.12" square 
boxes centered on each component, we found that the quality of the deconvolution 
improved as evidenced by comparison of the B-array maps to the A-array maps. 

When two point sources are barely resolved it is difficult to measure their flux 
densities independently. We attempted to use the AIPS task JMFIT to fit two gaus- 
sians to the A1-A2 double and found a strong covariance between the measured flux 
densities. Since we have already fitted for the flux densities of each component in 
the image deconvolution step, we instead measure the flux density of each image 
by summing the flux density of the CLEAN components in the box centered at 
the image position. This CLEAN flux density does not reliably measure the total 
flux density when the clean components do not reliably represent the total flux dis- 
tribution (e.g. for extended sources). Since all four images in MG0414 are nearly 
point-like, the CLEAN components contain essentially all of the flux and can be 
used directly to get the flux densities. Our flux density measurement was compared 
to the results of JMFIT by creating artificial data sets with the same Mv-sampling 
as typical data and four gaussian components with parameters measured from our 
A-array images. Gaussian white noise was added to all visibilities so that the noise 
in the image was comparable to the noise measured in the real images. The arti- 
ficial data sets were deconvolved in exactly the same way as the real data and the 
resulting images then measured with JMFIT and by our technique. The results 
are comparable for the case when all images are well resolved (A-array) but in the 
B-array data, the A1-A2 pair is marginally resolved and our flux densities for those 
images had a variance smaller than those of JMFIT by a factor of ~50. 

Weather conditions can make flux calibration at 15 GHz very difficult. The 
most important effects are varying absorption due to the proximity to the water 
absorption line at 22 GHz and wind-induced pointing errors of the individual array 
antennae. For example, we have measured variations in calibrator flux density of a 
factor of ~ 1.5 over a period of 25 minutes at the beginning of a snowstorm. (The 
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data from that particular observation are omitted from our analysis.) Since flux 
calibration errors can substantially affect the results of the analyses to be discussed 
below, these analyses are carried out using only those observations made when 
the wind speed was less than lOm/s, the object was above 30° elevation, there 
was no precipitation, and the phase stability of the atmosphere allowed reliable 
calibration of the complex antenna gains. Of the original 63 observations, three 
were lost to equipment failures at the VLA site, calibrations for two were lost to 
human error at operations, and one was lost to snow filling the dishes. Of the 
remaining 57 observations, we exclude from analysis five for high wind, eight for 
precipitation, two for atmospheric phase instability, and seven for low elevation, 
leaving 35 observations. Table || lists the measured flux densities for all components 
for each of the 35 observations. 
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Table 2: Table of measured flux densities (in mJy). 

Figure Q displays the measured light curves for all four images in MG0414. The 
variability is still clearly dominated by the flux calibration errors caused by variable 
atmospheric water content, differences in source and calibrator elevations, flux cali- 
brator modelling errors, and weather conditions. Even with these effects, the RMS 
variability in the resulting four light curves (without any adjustment for possible 
source variability) is only 3.3-3.9%. Since this includes both source variability and 
flux calibration errors, the relative flux calibration is at least this good over the 
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whole data set. 




Figure 4: Light curves for all components shown on a log scale so that shifting on 
the vertical axis is equivalent to scaling. 



4 Analysis 
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4.1 Variability in the Source 



Source variability is prerequisite to the determination of a time delay but the uncer- 
tainties in the flux calibration obscure any evidence for source variability in the light 
curves of Figure ||. In order to determine if there is source variability, we compare 
the variations in the ratios of various images. Figure || displays "ratio light curves" 
in which the flux density of each component is normalized to component B. For 
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Figure 5: Time history of the ratio of the flux density of each image to image B. 
Shown on a log scale so that shifting on the vertical axis is equivalent to scaling. 
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comparison, we generated an artificial data set with the same ui>-sampling as the 
real data, the average fluxes for the four components of MG0414 that we measure, 
and no source variability (created by by the AIPS task UVMOD). An estimate of 
the RMS noise (per visibility) was obtained from the phase calibrator observations 
in the real data and gaussian noise of this amplitude was added to the visibilities 
in the artificial data. This complete artificial data set was then processed with the 
same run-file as the real data. Table || displays the RMS variability in each flux 
ratio for both the real data and the artificial data described above. These results 
demonstrate that uncertainties in the map computation and flux estimation pro- 
cedures are comparable to the uncertainties in the calibration of the overall flux 
density scale. Thus, there is no evidence here that any source variability has been 
detected. 



Ratio Real Data Artificial Data 



A1:A2 0.9% 1.4% 

ALB 2.1% 1.9% 

A2:B 1.9% 1.4% 

A1:C 2.2% 2.6% 

A2:C 2.1% 2.8% 

B:C 2.9% 2.7% 



Table 3: Measured RMS variability in image flux density ratios for the data re- 
ported here and an artificial data set with similar instrumental properties. 



Another measure of variability is the first order structure function (Simonetti 
et al. 1985| ), 

V(r)= 1 -([s(t)-s(t + r)] 2 ) 

where s(t) is the image flux density at time t, r is the difference in time between 
two observations, and ( ) indicates the average.) We fit the structure function of the 
log of the radio data in decibels (referenced to 1 Jy), adopting a structure function 
of the form 

V(t) = Ct 1 - . 

We expect that the structure function of instrumental errors will be that of white 
noise, so we take the power law part of the structure function as an estimate of the 
maximum possible variability in the source and find 

, . . decibels, 

V(t) = (0.000139 — )r. 

day 

This implies that the signal-to-noise ratio for a single measurement of a change in 
flux density on time scales of interest (~20days) is at most ~0.3. 
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4.2 Time Delays and Correlated Variability 

The data presented above provide no evidence that we have detected variability in 
MG0414. However, it is possible that some of the variability in the data is due 
to the quasar and that such "real" variability would manifest itself as itself in the 
form of a correlation between flux density measurements separated by the lens time 
delay. To investigate this, we computed discrete correlation functions ( |Edelson~S 



Krolik 1988| ) for pairs of our light curves. No evidence for correlations was found in 



this analysis. 

A time delay detection technique with greater statistical efficiency than cross- 
correlation techniques has been developed by Press, Rybicki, and Hewitt ( |1992a ; 



1992b|) ; see also flRybicki fc Press 1992j ). Here we examine their x 2 statistic (here- 
after PRH-x 2 ). To apply their method, one assumes (or better, fits for) a flux 
ratio between two images and then, for each trial delay, creates a combined light 
curve that consists of the union of the two light curves with one of them shifted 
and scaled relative to the other. For each trial delay, we ask how likely (in a x 2 
sense) the combined light curve is to have originated from a process whose variabil- 
ity is described by the assumed structure function. We have performed Monte Carlo 
simulations by generating 500 trials using the sampling of the data from this work, 



the delays predicted by the Hewitt et al. (1992) models, the noise measured from 
the data (assumed to be white gaussian), and assuming the source variability to be 
generated by a gaussian process with the same structure function as that estimated 
from the data. For comparison, we also performed the same simulation assuming no 
gravitational lens (i.e. every image is an independent source of gaussian noise again 
with the same structure function as that estimated from the data but uncorrelated 
with the other images). Figure |6| displays a histogram of the absolute minimum of 
the PRH~x 2 curve on the interval [-45,45] for every image pair in the presence of 
the modelled gravitational lens. The simulations show that given the quality and 
quantity of the available data and the small (if present) variability of the source, 
extraction of a time delay is not possible by this technique. Only 22-31% of the 
measurements are within ±5d of the simulated delay. However, it is clear that the 
technique is detecting the presence of correlated variability since the correlated sim- 
ulations show 3-6 times more measurements within 5 d of the simulated delay than 
the uncorrelated (no gravitational lens) simulations. 

Even given the low signal-to-noise ratio of the present data, the PRH-x 2 method 
is interesting since it uses all available lag pairs and is thus a very sensitive way to 
measure correlated variability. Given sufficiently large data sets, one can expect to 
detect correlated variability even when the signal-to-noise ratio for a single lag pair 
is less than unity ( Hewitt et al. 1995| ) . 



From the measured MG0414 data, one can compute a distinct PRH-x 2 curve 
for every pair of images using the average flux ratio as an estimate of the true flux 
ratio. The six curves computed from the present data are displayed in Figure [7|. It is 
clear from examining the PRH-x 2 curves that there is a strong correlation between 
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FIGURE 6: Histograms of the absolute minimum of the PRH-x 2 ° n the interval 
[—45, 45] for 500 Monte Carlo trials assuming the presence of a gravitational lens. 
The model delay is marked with an arrow in each case. The bin on the left edge of 
the B vs. C plot is off-scale, its value is 121. 



the light curves at zero lag due to unmodelled errors in the absolute flux calibration 
and that none of them presents a convincing argument for any particular time delay. 
This is not surprising. The signal-to-noise ratio of 0.3 for the variability as measured 
in a single lag pair (see Section 4.1 above) can be compared to the case of Hewitt et 
gl ? s (li~995D MG1131+0456 data for which the same ratio is 0.7. Their simulations 
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of time delay observations in MG1131 show that 50 to 100 flux measurements are 
required for a measurement of the time delay. Clearly this study of MG0414 suffers 
because there are too few measured points, given the small signal-to-noise ratio. It 
is interesting to note that in the cases where the PRH-x 2 curves show a particular 
preference for the sign of the delay, it is consistent with the predictions of Kochanek's 
models 2, 3, and 4 (|Kochanek 1991) and the Hewitt et al. (1992) model. 



5 Conclusions 

Measurement of a time delay from these data is difficult at best because the variabil- 
ity is very small, the small size of the data set gives few lag pairs, and the zero lag 
correlation caused by errors in flux calibration dominates the PRH~x 2 curves. One 
could, in principle, modify the statistical model of the data used in the x 2 analysis to 
account for the zero-lag correlation of the flux density calibration errors. However, 
given the limited expectation of a determination of a time delay we postpone such 
analysis until a better data set is available. Our current analysis does not provide a 
convincing case for source variability or for any particular value of a time delay in 
MG0414. 

Improved experiments for measuring this time delay might use one or more of 
the following strategies: obtaining a larger data set, observing at a lower frequency 
where flux calibration is more precise, or (conversely) observing at a higher frequency 
where the source is more likely to be strongly variable. 
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Figure 7: PRH-x 2 curves for all pairs of components (69 degrees of freedom). Note 
the prominent feature at zero lag due to flux calibration errors. Also note that some 
of the curves, particularly those in the right column exhibit a strong preference for 
a particular sign of the delay. For ease of comparison, the curves are all plotted 
with the same vertical scale. 
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